begin
  
  model = "mcv"
  fi = addfile("../divvor_regional4_mcv.nc", "r")
;  fi = addfile("../power_divvor_"+model+"3.nc", "r")
  power_div  = fi->power_div
  power_vor  = fi->power_vor
  wavelength = fi->wavelength
  wavenumber = fi->wavenumber 

  data = new((/3, dimsizes(power_div)/), double)
  data(0,:) = power_div + power_vor
  data(1,:) = power_vor
  data(2,:) = power_div
  
  
  wks = gsn_open_wks("ps", "../power_div_vor_"+model+"4")
  res1                    = True
  res1@gsnFrame           = False
  res1@trXLog             = True
  res1@trYLog             = True
  res1@trXMinF            = 1
  res1@trXMaxF            = dimsizes(power_div)
  res1@trYMaxF            = 1.e+03
  res1@trYMinF            = 1.e-07
  res1@xyLineThicknesses  = (/4.0, 4.0, 4.0/)
  res1@tiXAxisString      = "Wavenumber (n)"
  res1@tiYAxisString      = "Kinetic Energy(k) (J kg~S~-1~N~)"
  res1@xyDashPatterns     = (/0, 0, 0/)
  res1@tmYROn             = True
  res1@tmYRBorderOn       = True
  res1@tmYRMajorOutwardLengthF = 0.0
  res1@tmYRMinorOutwardLengthF = 0.0
;  res1@tmBorderThicknessF = 1.0 
;  res1@tmBorderLineColor  = "grey"

  res2                = True
  res2@gsnFrame       = False
  res2@trXLog         = True
  res2@trYLog         = True
  res2@trXMinF        = min(wavelength(0:))
  res2@trXMaxF        = max(wavelength(0:))
;  res2@trYMinF        = res1@trYMinF
;  res2@trYMaxF        = res1@trYMaxF
  res2@trXReverse     = True
  res2@tmXTMode       = "Explicit"
  res2@tiXAxisString  = "Wavelength (km)"
  res2@tmXUseBottom   = False
  res2@tmXTLabelsOn   = True
  res2@xyDashPatterns = res1@xyDashPatterns
  res2@tmYROn   = res1@tmYROn
  res2@tmYRBorderOn = res1@tmYRBorderOn
  res2@tmYRMajorOutwardLengthF = res1@tmYRMajorOutwardLengthF
  res2@tmYRMinorOutwardLengthF = res1@tmYRMinorOutwardLengthF
;  res2@tmBorderThicknessF = res1@tmBorderThicknessF
;  res2@tmBorderLineColor = res1@tmBorderLineColor
  ;------------------------------------------------------------
  ; reference lines (power -3 and power -5/3)
  x1 = 4
  x2 = 4000
  y1 = 1.0e2
  y2 = exp(-3.0*log(x2/x1))*y1
  res1@tmXTOn       = False
  res1@xyLineColors = (/"grey"/)
  res1@xyLineThicknessF = 1.0
  res1@tiMainString   = ""
  plot = gsn_csm_xy(wks, (/x1,x2/), (/y1,y2/), res1)
  res1@tiMainString   = ""
  delete(res1@xyDashPatterns)
  res1@xyDashPatterns = 1
  y2 = exp(-5.0/3.0*log(x2/x1))*y1
  plot = gsn_csm_xy(wks, (/x1,x2/), (/y1,y2/), res1)
  ;------------------------------------------------------------
  ; kinetic energy spectrum
  delete(res1@xyLineColors)

  res1@xyLineColors       = (/"black", "blue", "red"/)
  res2@xyLineColors = res1@xyLineColors
  res1@xyDashPatterns = 0
;  plot = gsn_csm_x2y(wks, wavenumber2, wavelength2, E_lindborg, res1, res2)
  res1@tmXTOn         = False
;  plot = gsn_csm_xy(wks, wavenumber, power_div+power_vor, res1)
  print(dimsizes(ispan(11, dimsizes(wavenumber),1)))
  print(dimsizes(wavelength(10:)))
  plot = gsn_csm_x2y(wks, ispan(1, dimsizes(wavenumber), 1), wavelength(0:), data(:,0:), res1, res2)

  ;add legend
  lgres                     = True
  lgres@lgLabelFontHeightF  = 0.1
  lgres@vpWidthF            = 0.25
  lgres@vpHeightF           = 0.1
  lgres@lgPerimOn           = True
  lgres@lgItemOrder         = (/2,1,0/)
  lgres@lgDashIndexes       = (/0,0,0/)
  lgres@lgLineColors        = (/"black", "blue", "red"/)
  lgres@lgLineThicknessF    = 4.0
  labels                    = (/"total KE", "rotational KE", "divergent KE"/)
  legend                    = gsn_create_legend(wks, 3, labels, lgres)
  amres                     = True
  amres@amZone              = 0
  amres@amParallelPosF      = -0.25 ; left-right
  amres@amOrthogonalPosF    = 0.35 ; up-down
  annota = gsn_add_annotation(plot, legend, amres)
  
  draw(plot)
  frame(wks)
end

